package com.asa.numanaly;

import java.util.Arrays;

import com.asa.hutils.MathHUtils;
import com.asa.utils.Calculate;

/**
 * 线性方程组的直接求解
 * 
 * @author Administrator
 *
 */
public class F {

	/**
	 * 线性方程组求解 之前用过一个，直接调用了，懒！！！！
	 * 
	 * @param a
	 * @param b
	 * @return
	 */
	public static double[] fangchengzu1(double a[][], double b[]) {
		double[][] aremenber = remember_a(a);		// 因为二维数组的传递是地址传递，而这个方法又会改变原来地址的a的值
		// 所以我们所有的计算要和原有的a隔离开
		double[] bremenber = remember_b(b);
		double[] fangchengzudejie = Calculate.fangchengzudejie(aremenber, bremenber);
		return fangchengzudejie;
	}

	

	/**
	 * 部分选主元的Gauss高斯消去
	 * 这些东西是为了然截断误差比较小，所以需要选择主元，尽量让主元相互消去时的乘子比较大
	 * p322
	 * @param a
	 * @param b
	 * @return
	 */
	public static double[] fangchengzu2(double a[][], double b[]) {

		return null;
	}

	/**
	 * 具有比例因子部分选主元的Gauss高斯消去
	 * p324
	 * @param a
	 * @param b
	 * @return
	 */
	public static double[] fangchengzu3(double a[][], double b[]) {

		return null;
	}

	/**
	 * 矩阵的逆 这个是至今为止最好的方法，计算量不大
	 * 这个算法在书中没有写出来
	 * 但是我用例子把算法还原出来了
	 * p334
	 * @param a
	 * @param b
	 * @return
	 * @throws Exception
	 */
	public static double[][] fangchengzu000(double a[][]) throws Exception {



		double[][] result = new double[a.length][a.length];
		for (int i = 0; i < a.length; i++) {
			double[] b = new double[a.length];
			b[i] = 1;
			result[i] = fangchengzu1(a, b);


		}

		return result;
	}

	private static double[][] remember_a(double[][] a) {
		double[][] aremenber = new double[a.length][a.length];
		for (int i = 0; i < aremenber.length; i++) {
			for (int j = 0; j < aremenber.length; j++) {
				aremenber[i][j] = a[i][j];
			}
		}
		return aremenber;
	}

	private static double[] remember_b(double[] b) {
		// TODO Auto-generated method stub
		double[] bremenber = new double[b.length];
		for (int i = 0; i < b.length; i++) {
			bremenber[i] = b[i];
		}
		
		
		
		return bremenber;
	}
	
	
	
	
	/**
	 * LU分解
	 * 书上的算法写得有点感人
	 * 我找了个，核心原理是一样的
	 * p350
	 * @param a
	 * @return
	 * @throws Exception
	 */
	public static double[][][] fangchengzu4(double a[][]) {

//		a = new double[][]{{1,1,0,3},{2,1,-1,1},{3,-1,-1,2},{-1,2,3,-1}};

		double[][][] result = new double[2][][]; 
		
		
		double u[][] = new double[a.length][a.length];
		double l[][] = new double[a.length][a.length];
		for (int k = 0; k < a.length; k++) {
			
			l[k][k] = 1;//怀疑这行

			for (int j = k; j < a.length; j++) {
				double asa = 0;
				for (int s = 0; s <= k-1; s++) {
					asa += l[k][s]*u[s][j];
				}
				u[k][j] = a[k][j] - asa;
				
				if (k==j) {
					System.out.println(k + " " + j);
				}
				
			}
			
			
			for (int i = k+1; i < l.length; i++) {
				double all = 0;
				for (int s = 0; s < k; s++) {
					all+=l[i][s]*u[s][k];
				}
				l[i][k] = (a[i][k]-all)/u[k][k];
			}
			
		}
		
		for (int i = 0; i < l.length; i++) {
			for (int j = 0; j < l[i].length; j++) {
				System.out.print("\t"+l[i][j]);
			}
			System.out.println();
		}
		System.out.println(".......................");
		for (int i = 0; i < u.length; i++) {
			for (int j = 0; j < u[i].length; j++) {
				System.out.print("\t"+u[i][j]);
			}
			System.out.println();
		}

		result[0] = l;
		result[1] = u;

		
		return result;
	}
	
	
	/**
	 * LDLt分解
	 * p360	也是找的其他资料
	 * Kincaid_Cheney_Numerical_Mathematics_and_Computing__Sixth_Edition 	p327
	 * 
	 * @param a
	 * @return	中间的那个矩阵是对角矩阵，对角线以外全为0，对角线上的是特征值(对角线上是特征存疑)
	 */
	public static double[][][] fangchengzu5(double a[][]) {
		double[][][] result = new double[2][][]; 

		double[][] l = new double[a.length][a.length];
		double[] d = new double[a.length];
		d[0] = a[0][0];
		
		
		for (int j = 0; j < d.length; j++) {
			
			l[j][j] =1;
			
			double asa = 0;
			for (int v = 0; v < j; v++) {
				asa += d[v]*l[j][v]*l[j][v];
			}
			System.out.println(asa+"..........");
			d[j] = a[j][j] -asa;
			System.out.println(d[j]);

			for (int i = j; i < d.length; i++) {
				l[j][i] = 0;
				asa = 0;
				for (int v = 0; v < j; v++) {
					asa += l[i][v]*d[v]*l[j][v];
				}
				
				
				l[i][j] = (a[i][j] -asa)/d[j];
				
				
			}
			
		}
		
		
		for (int i = 0; i < l.length; i++) {
		for (int j = 0; j < l[i].length; j++) {
			System.out.print("\t"+l[i][j]);
		}
		System.out.println();
	}
		
		
		for (int i = 0; i < d.length; i++) {
		System.out.println(d[i]);

	}
		
		
		result[0] = l;
		double[][] D = new double[l.length][l.length];
		for (int i = 0; i < d.length; i++) {
			 D[i][i] = d[i];
		}
		result[1] = D;

		
		return result;
		
		
	}

	
	/**
	 * 	特征值分解
	 * 	那本书上面没有，我补充的
	 * 	这里只返回 特征矩阵X 和特征值组成的对角矩阵A
	 * 	特征矩阵的逆，就自己再去求吧
	 * 
	 * @param a 三维的向量，其中0是特征向量组成的矩阵X，1是特征值组成的矩阵A
	 * @return
	 * @throws Exception 
	 */
	public static double[][][] XAX(double a[][]) throws Exception {
		double[] juzhengtezhengzhi = MathHUtils.juzhengtezhengzhi(a);
//		Arrays.sort(juzhengtezhengzhi);
		double[][] juzhengtezhengxiangliang = MathHUtils.juzhengtezhengxiangliang(a,juzhengtezhengzhi);
//		for (int i = 0; i < juzhengtezhengxiangliang.length; i++) {
//			for (int j = 0; j < juzhengtezhengxiangliang[i].length; j++) {
//				System.out.print(juzhengtezhengxiangliang[i][j]+"   ");
//			}
//			System.out.println();
//		}
		
		
//		for (int i = 0; i < juzhengtezhengzhi.length; i++) {
//			System.out.print(juzhengtezhengzhi[i]+"   ");
//		}
		
		double[][] A = new double[juzhengtezhengzhi.length][juzhengtezhengzhi.length];
		
		for (int i = 0; i < A.length; i++) {
			for (int j = 0; j < A[i].length; j++) {
				if (i!=j) {
					A[i][j] = 0;
				}else {
					A[i][j] = juzhengtezhengzhi[i];
				}
			}
		}
		
		
		double[][][] result = new double[2][][];
		result[0] = juzhengtezhengxiangliang;
		result[1] = A;
		
		return result;
	}
	
	
	
	
	
	/**
	 * SVD分解	奇异值分解
	 * 	那本书上面没有，我补充的
	 * 
	 * @param a
	 * @return	返回的是长度为三，其中第一个是S，第二个是特征，第三个D
	 * @throws Exception 
	 */
	public static double[][][] SVD(double a[][]) throws Exception {
		
		double[][] zhuanzhi = MathHUtils.zhuanzhi(a);
		
		double[][] chenfa1 = MathHUtils.chenfa(a, zhuanzhi);
//		for (int i = 0; i < chenfa1.length; i++) {
//			for (int j = 0; j < chenfa1[i].length; j++) {
//				System.out.print(chenfa1[i][j]+"   ");
//			}
//			System.out.println();
//		}
		
		double[][] chenfa2 = MathHUtils.chenfa(zhuanzhi,a);
		
//		for (int i = 0; i < chenfa2.length; i++) {
//			for (int j = 0; j < chenfa2[i].length; j++) {
//				System.out.print(chenfa2[i][j]+"   ");
//			}
//			System.out.println();
//		}
		
		
		double[][] juzhengtezhengxiangliang = MathHUtils.juzhengtezhengxiangliang(chenfa1);
		
		
//		for (int i = 0; i < juzhengtezhengxiangliang.length; i++) {
//			for (int j = 0; j < juzhengtezhengxiangliang[i].length; j++) {
//				System.out.print(juzhengtezhengxiangliang[i][j]+"   ");
//			}
//			System.out.println();
//		}
		double[] juzhengtezhengzhi = MathHUtils.juzhengtezhengzhi(chenfa1);
		
//		for (int i = 0; i < juzhengtezhengzhi.length; i++) {
//			System.out.println(juzhengtezhengzhi[i]);
//		}
		
		
		double[][] juzhengtezhengxiangliang2 = MathHUtils.juzhengtezhengxiangliang(chenfa2);
//		for (int i = 0; i < juzhengtezhengxiangliang2.length; i++) {
//			for (int j = 0; j < juzhengtezhengxiangliang2[i].length; j++) {
//				System.out.print(juzhengtezhengxiangliang2[i][j]+"   ");
//			}
//			System.out.println();
//		}
		
		
		double[][] A = new double[juzhengtezhengzhi.length][juzhengtezhengzhi.length];
		
		for (int i = 0; i < A.length; i++) {
			for (int j = 0; j < A[i].length; j++) {
				if (i!=j) {
					A[i][j] = 0;
				}else {
					A[i][j] = Math.pow(juzhengtezhengzhi[i], 0.5);
//					System.out.println(A[i][j]);
				}
			}
		}
		
		double[][][] result = new double[3][][];
		
		result[0] = juzhengtezhengxiangliang;
		result[1] = A;
		result[2] = juzhengtezhengxiangliang2;

		return result;
	}

	
	
	
	/**
	 * LLt分解
	 * p361	也是找的其他资料
	 * 
	 * @param a
	 * @return
	 */
	public static double[][] fangchengzu6(double a[][]) {

		double[][] l = new double[a.length][a.length];
		
		for (int k = 0; k < l.length; k++) {
			
			double asa = 0 ;
			
			for (int s = 0; s < k; s++) {
				asa = asa + l[k][s]*l[k][s];
			}
			
			
			
			l[k][k] = Math.pow(a[k][k] - asa,0.5);
			for (int i = k+1; i < l.length; i++) {
				
				asa = 0;
				for (int s = 0; s < k; s++) {
					asa = asa + l[i][s]*l[k][s];

				}
				
				l[i][k] = (a[i][k] - asa)/l[k][k];
				
			}
			
		}
		for (int i = 0; i < l.length; i++) {
			for (int j = 0; j < l[i].length; j++) {
				System.out.print("\t"+l[i][j]);
			}
			System.out.println();
		}
		
		return l;
		
		
		
	}
	
	
	/**
	 * 对三角线性方程组的Crout分解
	 * p364	也是找的其他资料
	 * 
	 * @param a	注意这是增广矩阵注意这是增广矩阵注意这是增广矩阵注意这是增广矩阵
	 * @return
	 */
	public static double[] fangchengzu7(double a[][]) {
		
		
		
		
		double u[][] = new double[a.length][a.length];
		double l[][] = new double[a.length][a.length];
		double z[] = new double[a.length];
		double x[] = new double[a.length];

		for (int i = 0; i < x.length; i++) {
			u[i][i]=1;
		}
		
		
		l[0][0] = a[0][0];
		u[0][1] = a[0][1]/l[0][0];
		z[0] = a[0][a.length]/l[0][0];
		
		for (int i = 1; i < a.length-1; i++) {
			
			l[i][i-1] = a[i][i-1];
			l[i][i] = a[i][i] -l[i][i-1]*u[i-1][i];
			u[i][i+1] = a[i][i+1]/l[i][i];
			z[i] = (a[i][a.length]-l[i][i-1]*z[i-1])/l[i][i];
			
		}
		
		l[a.length-1][a.length-2] = a[a.length-1][a.length-2];
		l[a.length-1][a.length-1] = a[a.length-1][a.length-1] - l[a.length-1][a.length-2]*u[a.length-2][a.length-1];
		z[a.length-1] = (a[a.length-1][a.length] - l[a.length-1][a.length-2]*z[a.length-2])/l[a.length-1][a.length-1];
		
		
		x[a.length-1] = z[a.length-1];
		for (int i = 0; i < z.length; i++) {
			System.out.println(z[i]);
		}
		System.out.println(".......................");

		for (int i = x.length-2; i >=0; i--) {
			x[i] = z[i] - u[i][i+1]*x[i+1];
		}
		System.out.println(".......................");

		for (int i = 0; i < x.length; i++) {
			System.out.print(x[i]+"   ");
		}
		System.out.println(".......................");

		for (int i = 0; i < l.length; i++) {
			for (int j = 0; j < l[i].length; j++) {
				System.out.print("\t"+l[i][j]);
			}
			System.out.println();
		}
		System.out.println(".......................");
		for (int i = 0; i < u.length; i++) {
			for (int j = 0; j < u[i].length; j++) {
				System.out.print("\t"+u[i][j]);
			}
			System.out.println();
		}
		return x;
		
	}
	
	
	
	
	
	
	
	
	
	public static void main(String[] args) throws Exception {
//		double a[][] = { { 1, 2, -1 }, { 2, 1, 0 }, { -1, 1, 2 } };
//
//		double[][] juzhengni = fangchengzu000(a);
//		juzhengni = MathHUtils.zhuanzhi(juzhengni);
//
//		for (int i = 0; i < juzhengni.length; i++) {
//			for (int j = 0; j < juzhengni.length; j++) {
//				System.out.print(juzhengni[i][j] + "    ");
//			}
//			System.out.println();
//		}
//		System.out.println("==========================");
//		for (int i = 0; i < a.length; i++) {
//			for (int j = 0; j < a.length; j++) {
//				System.out.print(a[i][j] + "    ");
//			}
//			System.out.println();
//		}
//		double a[][] = {{1,1,0,3},{2,1,-1,1},{3,-1,-1,2},{-1,2,3,-1}};
//		double a[][] = {{4,-1,1},{-1,4.25,2.75},{1,2.75,3.5}};
//		fangchengzu4(a);
//		double a[][] = {{4,-1,1},{-1,4.25,2.75},{1,2.75,3.5}};
//		double[][][] fangchengzu5 = fangchengzu5(a);
//		double[][] L = fangchengzu5[0];
//		double[][] D = fangchengzu5[1];
//		for (int i = 0; i < L.length; i++) {
//			for (int j = 0; j < L[i].length; j++) {
//				System.out.print(L[i][j] + "    ");
//			}
//			System.out.println();
//		}
//		for (int i = 0; i < D.length; i++) {
//			for (int j = 0; j < D[i].length; j++) {
//				System.out.print(D[i][j] + "    ");
//			}
//			System.out.println();
//		}
		
		
		
		
		
		
//		fangchengzu6(a);
		
		
//		double a[][] = {{2,-1,0,0,1},{-1,2,-1,0,0},{0,-1,2,-1,0},{0,0,-1,2,1}};
//		
//		fangchengzu7(a);
		
		
		
		
		
//		double[][] a = {{-1,0,0,1,0},{0,-1,0,-1,1},{0,0,-1,1.0/2,-1},{0,0,0,-1.0/6,0.5},{0,0,0,1.0/24,-1.0/6}};
//		double[] b = {1,-0.5,1.0/6,-1.0/24,1.0/120};
//		double[] fangchengzu1 = fangchengzu1(a, b);
//		
//		for (int i = 0; i < fangchengzu1.length; i++) {
//			System.out.println(fangchengzu1[i]);
//		}
		
		
//		double x[] = new double[b.length];
//		double sum = 0;
//		double xmult = 0;
//		for (int k = 0; k < b.length-1; k++) {
//			for (int i = k+1; i < b.length; i++) {
//				xmult = a[i][k]/a[k][k];
//				a[i][k]=xmult;
//				for (int j = k+1; j < b.length; j++) {
//					a[i][j]=a[i][j]-(xmult)*a[k][j];
//				}
//				b[i]=b[i]-(xmult)*b[k];
//			}
//		}
//		x[b.length-1] = b[b.length-1]/a[b.length-1][b.length-1];
//		for (int i = b.length-2; i >=0; i--) {
//			sum = b[i];
//			for (int j = i+1; j < x.length; j++) {
//				sum = sum - a[i][j]*x[j];
//			}
//			x[i] = sum/a[i][i];
//		}
//		
//		
//		
//		System.out.println("解：");
//		for (int i = 0; i < a.length; i++) {
//			System.out.println(x[i]);
//		}
		
//		System.out.println("============================");
//		for (int i = 0; i < b.length; i++) {
//			System.out.println(b[i]);
//		}
//		System.out.println("============================");
//		for (int i = 0; i < a.length; i++) {
//			for (int j = 0; j < a.length; j++) {
//				System.out.print(a[i][j]+"\t");
//			}
//			System.out.println();
//		}
		
		
		
		
		
		
		//SVD
//		https://www.cnblogs.com/endlesscoding/p/10033527.html
//		double[][] a = {{1,5,7,6,1},
//				{2,1,10,4,4},
//				{3,6,7,5,2}};
//
//		SVD(a);
		
		
		
		System.out.println("========");
		
		
		double[][] a = {{2,1,1},
		{1,3,2,},
		{1,2,2}};
		
		
		double[][][] doolittle_LR = Doolittle_LR(a);
		
		
		
	}
	
	
	
	
	
	/**
	 * 使用Doolittle方法求得LR分解
	 * @param a 要分解的矩阵
	 * @return	0位置是L，1位置是R
	 */
	public static double[][][] Doolittle_LR(double[][] a) {
		
		

		
		double[][] l = new double[a.length][a.length];
		double[][] r = new double[a.length][a.length];
		
		
		for (int i = 0; i < l.length; i++) {
			l[i][i] = 1;
		}
		
		
//		for (int i = 0; i < r.length; i++) {
//			for (int j = i; j < r[i].length; j++) {
//
//				double asa = 0;
//				for (int k = 0; k < i-1; k++) {
//					asa += l[i][k] * r[k][j];
//				}
//				r[i][j] = a[i][j] - asa;
//				
//			}
//		}
		
		
		
		for (int k = 0; k < r.length; k++) {
			
			for (int i = k; i < r.length; i++) {
				double asa = 0;
				
				for (int p = 0; p < k; p++) {
					asa += l[k][p] * r[p][i];
				}
				r[k][i] = a[k][i] - asa;
	
			}
			
			
			for (int i = k+1; i < r.length; i++) {
				double asa = 0;
				for (int p = 0; p < r.length; p++) {
					asa += l[i][p] * r[p][k];
				}
				l[i][k] = (a[i][k] - asa)/r[k][k];
			}
			
			
			
		}
		
		
		
		
		
		System.out.println("L=");
		MathHUtils.print(l);
		System.out.println("R=");
		MathHUtils.print(r);
		
		
		double[][][] result = new double[2][][];
		result[0] = l;
		result[1] = r;
		
		
		return result;
		
		
		
	}
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	

}
